#!/bin/bash
#
# Script to analyze the MSD, must be run after the calculation.

# Extract info from log file
awk /"^   Step "/,/"^Loop time"/ log.lammps|tail -10003 > Info.dat

# Time step and offset info.
dt=0.001 # ps
s0=`head -2 "Info.dat"|tail -1|awk '{print $1}'`
m0=`head -2 "Info.dat"|tail -1|awk '{print $10}'`

# Generate the gnuplot script for visualization
cat > .plot.script << EOF
set term png
set out 'msd.png'

set xlabel 't (ps)'
set ylabel 'msd (Angstrom^2)'

#set yr [0:*]
unset key

f(x) = a * x + b
fit f(x) 'Info.dat' u ((\$1-${s0})*${dt}):(\$10-${m0}) via a, b

lbl_text = sprintf("D = %.4g Ang^2/ps", a/6)
set label lbl_text at 100,2500

plot 'Info.dat' u ((\$1-${s0})*${dt}):(\$10-${m0}) w p, f(x) w l lw 2

set term post enha color 20
set out 'msd.eps'
rep
EOF
# Draw the curve
gnuplot .plot.script
evince msd.eps

# Cleanup
#rm -rf .plot.script fit.log #Info.dat
